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Abstract 

In this paper, second order Godunov methods are reviewed. The early versions by Colella 
and Woodward (PPM) and van Leer (MUSCL) are described in their original form. The 
simplification of these by Roe, based on approximate Riemann solver, is then presented. 
Attention is next given to the improvement in MUSCL due to Hancock and van Leer leading 
to a fuller paper by Huynh. Finally, brief reference is made to TVD and ENO schemes due 
to Harten. 
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1. INTRODUCTION 


Godunov’s method for solving unsteady problems in Gas Dynamics is described and dis- 
cussed at length in Holt (1984). In the succeeding decade several extensions of the method have 
been proposed which increase its accuracy to second order while retaining the properties in the 
original method of Monotonicity and absence of oscillations in shock capturing. The contribu- 
tions most clearly related to the first Godunov scheme, notably by Colella and Woodward 
(1984), Roe (1986) and van Leer (1979) are reviewed and added to the earlier account. At the 
time of writing Holt (1984), these contributions had either not been completed or not yet recog- 
nized as Godunov extensions. The revised chapter in my book (to appear as a 3rd edition) will 
include these extensions and the present report is a preliminary version of this coverage. 

The extensions treated specifically are the MUSCL (Monotonic Upstream-centered Scheme 
for Conservation Laws) scheme of van Leer, the Piecewise Parabolic Method (PPM) of Colella 
and Woodward and the characteristic based scheme of Roe. This report firstly deals with these 
methods as originally presented. Thus PPM is described only in application to the unidirectional 
wave equation, with indication of its extension to the Eulerian Gas Dynamic Equations. The 
MUSCL scheme, as developed by van Leer (1979) is applied to the one dimensional Lagrangian 
equations. Roe’s scheme, freer than the other schemes of algebraic details, is presented is gen- 
eral form. It is hoped that in the final version of the Godunov chapter both Colella and van Leer 
will provide me with versions of their respective methods which are easier for Graduate Students 
to understand. This final version will (as it should) cover TVD and ENO schemes. 

1. Godunov Extensions The method proposed by Godunov (1959) for solving problems of one 
dimensional gas Dynamics is of first order accuracy and has the important property of monotoni- 
city. This requires that when the method is applied to an initial value problem in which the unk- 
nown has monotonic behavior at the outset, then the calculated values of the unknown at all later 
times remain monotonic in character. The monotonic property is crucial when dealing with 
compression waves, to ensure that the only shock waves which develop result from local 
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collapse of continuous compressions and exclude spurious socks resulting from faults in the 
numerical scheme. 

In its original form the method was applied to one dimensional shock tube problems and 
generalized to calculate steady flow past two and three dimensional bodies, including one of the 
earliest solutions to the blunt body problem. These applications provided valuable information 
for basic engineering analysis but the first order accuracy limited detailed flow field analysis. 

Godunov proved that his original scheme was restricted to first order accuracy and could 
not be extended to higher order without sacrificing the monotonicity property. Godunov (with 
colleagues, see Alalykin et al (1970)) then proposed a second order scheme, of predictor- 
corrector type, applied over successive half time steps. The predictor uses an implicit formula 
extending over three space steps and is carried out by a double sweep process extending between 
two outside boundaries. The data from the predictor step are smoothed before being used in the 
corrector step, which is analogous to Godunov’s first order scheme. Monotonicity is thus 
preserved, but applications have been limited to one dimensional problems and any shock waves 
or shock type discontinuities are fitted so that the scheme is tied to a moving network. Further 
research on this scheme is recommended. 

In applying Godunov’s first scheme to unsteady one dimensional gas flow, the Eulerian 
equations of motion are used. At a given time the disturbed flow regime is divided into cells 
(usually of equal width) and the unknowns, determined from the scheme up to that time, are 
represented as staircase formations with constant values in each cell (averages of current values 
at the left and right boundaries of the cell). The solution is advanced by solving successive 
Riemann problems in adjacent cells. The staircase representation is responsible for the limita- 
tion of Godunov’s first method to first order, although it also automatically carries monotonic 


character. 
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Several authors attracted by Godunov’s first method, because of its relevance to the physics 
of problems in gas dynamics, have proposed improvement in its accuracy by removing the res- 
triction that mean values of unknowns in cells should be constant. Instead they propose that 
unknowns should be represented by linear or higher order functions of distance across the cell, 
still solving Riemann problems at all boundaries. Such a feature is common to all so-called 
higher order Godunov methods identified principally with Colella, Roe and van Leer. 

Colella’s approach is an outgrowth of the Piecewise Parabolic Method (PPM) applied to the 
one directional wave equation (Colella & Woodward (1984)). Roe (1986) uses a generalized 
solution of Riemann problems for linear conservation laws while van Leer bases his treatment on 
improvements of his MUSCL scheme (van Leer (1979)). 

We shall describe all these extensions of Godunov methods and also refer to two devices 
for controlling monotonicity, the Total Variation Diminishing and ENO Schemes of Harten 
(Davis (1984)). 

2. Piecewise Parabolic Method. 


Following Colella & Woodward (1984) we consider the one directional wave equation 


du du A 
— + a— = 0 
3t 


u(£, 0) = u 0 (£) 


( 2 . 1 ) 


This equation is discretized - with equal space and time step Ai;, and At. The disturbed distribu- 
tion of u is divided into equal cells j, j + 1 -- with intercell boundaries) + V 2 , as shown 



£j + y, is the boundary between cells jand jt-1 
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Knowing u £ the value of u (£j ) at time t n , we wish to calculate u" + 1 The average value of u 
between £j + y 2 and E,j_ \ h at t = t n is 

<2 - 2> 

^ ^j — ^j+ Vi ~ ^ j- x h 

Stability requires a At < A £ in difference calculations. Following van Leer (1979) we construct a 
piecewise polynomial interpolation for u (£,) satisfying 



1 f 

■* 5 j-% 


u(£) 


To apply this at time t = t n + At we substitute the exact solutions of Eq. (2.1) 


u!' +1 = _Lj^»»„(4-aAt) d 5 
3 Ac, j J s 

and represent the integrand uniquely by a parabolic profile 


(2.3) 


“(£) = u L,j+ x 


AUj+u 6 j(l -x) 


(2.4) 


5j- . /2 < ^ < ^ j+ 1/2 

The coefficients in this quadratic expression are uniquely determined in terms of u" and limiting 
values at the left and right sides of boundary ^ + Vz 



lim u (£) = u L j 

£ (4d- w)+ 


lim u (£) = u R . 

% —■ ► (4 j+ &)- 


(2.5) 
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Au j= u R,j~ u L.j> u 6,j = 6 


“j - -Z ( u L,j + u R,j> 


To calculate u L and u R z we use an approximation for Uj + \ h which does not lie outside the range 
(Uj, Uj + ,). At smooth parts of the solution, away from extreme 


U L. j+ 1 _ U R, j“ u j+ Vi 

At all interfaces, if necessary, u L j and u R j are modified to ensure that u remains monotonic in 
(^j_ i h , , /2 ). Given averages of u " of u in nearby zones, we interpolate a value of u i/z from 

the indefinite integral of u, calculated at five consecutive points, fitted by a quartic which is dif- 
ferentiated at ^ j + 1 / 2 . We then find 

Uj+ 1/2 = 12 + U ^ + 1} " ~12 (u * 2 + l} (2 ' 6) 
In general u L j+1 and u R 0 -are equal to u j+ , /2 . If the interpolation function gives values of u j+ ^ 
outside the range (u L j+1 u R J then the values of the latter pair need to be reset to satisfy mono- 
tonic behavior; details are given in Colella and Woodward (1984). 

In the final step to get a formula for u " + 1 we introduce averages of the interpolation func- 
tions 


f U. L(y 


f j* ’/!, R<y)- y J&;*'*©* 


These evaluate to 


f " + I/s. L(y) = U R. j - J j Au j - (1 - b x) u 6ij \ for x = y/A^j 


(2.7) 


(2.8) 



(2.9) 


f n 

T j+ '/2, 


R(y) = u L, 



Au 


j+i 



for x = y/ 


and finally 


Uj +1 = u" + a 


At 


u 


y x h ^ji- l /2 


where 


u j+ 1 / 2 = f j+ Vi, L^) if a>0 
= f£.*/ 2 , r (-aAt) if a < 0 . 

Generalizations of this scheme, applied to the one dimensional Lagrangian and Eulerian 
equations include the following steps 

1 . Interpolation of the initial dependent variable distribution. 

2. Use of characteristic equations to find dependent variable values on each side of zone 
boundaries. 

3. Solution of Riemann problems at boundaries to determine numerical fluxes there. 

3. The MUSCL Scheme 

Van Leer (1979) introduced a higher order Godunov method which preceded the PPM 
method and also led to his flux or fluctuation splitting schemes. Its original version, called 
Monotonic Upstream-centered Scheme for Conservation Laws (MUSCL), was written by Wood- 
ward and later refined in collaboration with van Leer. 

The MUSCL scheme analysis is described by van Leer in terms of the Lagrangian form of 
the one dimensional unsteady flow equations to be mapped subsequently onto an Euler grid. 
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van Leer writes the Lagrangian equations 


3V (^“u) _ 

at " a^ 


o 


^ + I »|t = F 


at 




(3.1) 


(3.2) 


dE d(x«u £l = uF + G 

at a^ 


ax 

at 


= u 


(3.3) 


(3.4) 


Here t is time, £, is the mass coordinate related to the space coordinate x and volume coordinate 

Xby 


d& = V _1 x a dx = V' 1 ^- = V" 1 dx (3.5) 

a+ 1 

0 plane 

a = 1 cylindrical symmetry 
2 spherical 

V, u, E and p denote specific volume, velocity, total energy and pressure respectively. F and G 
represent source terms. If the specific internal energy is e 


E = e + 


2 


The equation of state is 


P = P(V, e) 


(3.6) 


(3.7) 


For an ideal perfect gas this simplifies to 


P = (7“ 


(3.8) 
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The Lagrangian sound speed is C where 


C 2 = - 


iL 

av 


adiabatic 


v 


related to the spatial sound speed by 


c =v 


The above equations can be written in characteristic form 


(3.9) 


(3.10) 


, T _ , „_i . i auVC 

dJ = du - C dp =■ 


+ F - C -, G 


de 


•dt 


on — = -x a C, r~ characteristic 
dt 


(3.11) 


dJ + = du + C _i dp=^-^^+F + CTG 


9p 

de 


-dt 


(3.12) 


on 


— = x a C, r 4- characteristic 


dt 


dS = dp + C 2 dV = G 


d£ 

de 


dt 


(3.13) 


on — - 0, streamline 
dt 

Across a shock wave the following jump equations are derived from the integral form of the 
above governing equations. 


±W(V*~ V) + (u*-u) = 0 


(3-14) 


± W (u* - u) - (p* - p) = 0 


(3.15) 
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± W (E* - E) - (u*p* - up) = 0 


( 3 . 16 ) 


Post shock values are indicated by an asterisk. The Lagrangian shock speed is ± W. 


Discretization 


The disturbed region is divided into cells of thicknesses A£. Interface values are denoted 
by integers. Values at a cell center and average values there (barred) are denoted by half 
integers. The notation is shown in Table 1. 


Svmbol 


& + ^ 

Xi + 1 

A, + i/2(t A i + 1/2 X 

t l 


Qi+ 1 / 2 1 Q l+1 ^ 2 
Qi+1/2, Q' +1/2 
Qi,Q l 
<Q>i 
< Q>i+ 1/2 
AiQ,A'Q 
Ai+l/2<3/ A, + i/ 2 ( 
Aj+i/2<3/A i+1/ f 2 y 


TABLE 1 

Notation Used in the Grid 
Definition 

Mass, Euler, volume coordinate of zone boundary " ~~ 

|(£ t + (»+i )■ mass-averaged mass coordinate of zone (&,&+ 1 ) 

2 (Xi + Xt+i)i volume-averaged volume coordinate of zone (x»>X;+i) 

(i + l (i , Vt + 1 Xi 

Initial time level 

t° + At, final time level 

Mass-averaged value of Q in zone ((i,(i+i) at t°, t 1 
Volume-averaged value of Q in zone (x,-,x i+ i) at t 0 ,/ 1 
Value of Q at the boundary £, at t 0 ,! 1 
Average value of Q at the boundary during time step 
Average value of Q in zone (&,&+,) during time step 

Qi+1/2 - Qi-1/2, Q i+1/2 - Q '~ 1/2 

Mass-averaged value of dQ/d^ in zone (&,&+ 1 ) at t° 

Volume-averaged value of dQ/d\ in zone (*;, y,+i) at t° 


We extract 


^\+Vi ^ + 1 + 'h ^ ^i+1 

In each slab, at the initial time t° we approximate initial values of V, u and E by linear distribu- 
tions e.g. 


V(tU) = V i + 1/2 + 


A i + ./:V 


A i + Vi ^ 




( 3 . 17 ) 


where 
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5m 

v i+ ./ 2 = (A i+ ./ 2 ^r 1 J v(t 0 ,^ 


Aj + ./ 2 V 

A i + '/ 2 £ 



= A i + I/2 V(t°, 

We retain only V, u, E, AV, Au and AE and these are to be updated. 
We can replace E by p. 

In practice we use 


( 3 . 18 ) 


( 3 . 19 ) 


— 1 2 

e i + '/2 = Ej + 1/2 - ~Uj + i/ 2 + 0 


m 2 


P, + >/ 2 = p(Vi + vs , ej + vs) + o ^ (A£) 2 


( 3 . 20 ) 


Note that slope values are independent of slab averages. Updated slab averages are derived by 
integrating the basic equations. 


V l * w = V i+1/4 + 


At 


A i+'/£ 


(<x a U> i+1 - <x a u> i ). 


( 3 . 21 ) 


u ,+,/2 = U|+vs “ T^T ( <X V i+l ~ < x Vi) 


A i+'/ 2 ^ 


( 3 . 22 ) 


+ (<apV/x> i+1/2 + <F> i+l/2 ) At, 


E 1+ ' /2 = Ej + vs - i At t (<x a up> i+1 - <x“up >j) + 


A i+'/2^ 


( 3 . 23 ) 


(<uP> i+1/! + <G>, + , /2 ) At; 


These hold even when discontinuities exist within the slab. Time averages appearing in these 
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equations need to be estimated with first order accuracy. 

To update AV, Au and Ap we need to estimate interface values of V, u and p at the end of 
this time interval and take their differences. For example 

A i+1/2 V = V i+I - V (3.24) 

Limiting values at interfaces 
These are given by 

Vi± = V i±K ±y4 i± *V 

u i±=u j±1 4± yA i±l/2 u (3.25) 

, 1 T 

Pi± ~ Pi ±'A - Ai + '/2 P 

i+ right side i- left side. 

With these discontinuous values defined at interfaces we solve breakdown formulae (simplified 
from those used in Godunov’s first scheme) to find new values at interfaces. These are denoted 
by u * , p * , etc. 

We also calculate the time derivatives. 


du 

* 

V 

* 

and 

av 

at 

W J 

i 

at 

i 

at ^ 


from the equations of motion in characteristic form. 


Details of Lagrangian scheme. 


( 1 ). Compute V ± , u ± and p+ — » C+ etc. 
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(2) . Find u* , p* , W+ V+. 

(3) . get (3u/dt)* (dp/3t)* (dV/dt)* from characteristic eqs. 

(4) . Integrate through the half time step w.r.t. time. Then interpolate. 

(5) . Use full time step with averaged value. 

(6) . Remap on to the Eulerian grid. 

Full formulae are given on pp. 1 12 and 1 13 in van Leer (1979). 

4 . An improved MUSCL scheme for the Euler equations. 

The principal obstacle to the extension of Godunov’s original method for the Euler equa- 
tions to second order accuracy is the preservation of monotonicity or monotonic property, in the 
difference scheme used. This is accentuated when such extensions are based on shock capturing 
rather than shock fitting, a feature common to virtually all the higher order schemes. 

To avoid breakdown of monotonic property a constraint must be introduced when, at any 
stage in the computation, the distribution in a calculated quantity indicates a local discontinuity 
or local extremum. Much effort has been spent on the identification of the optimum constraint 
which will, at the same time, ensure monotonic properties and maintenance of second order 
accuracy. 

A recent paper by Huynh (1995) reviews the history of constraint improvement and pro- 
poses a new constraint which appears to be closer to the ideal required. He incorporates this in a 
new approach to the numerical solution of the one dimensional Euler equations by a second 
order Godunov method. A basic feature of this is the simplification of the second order upwind 
scheme in MUSCL due to Hancock (1980). We begin, with his description of this for the simple 
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mass convection equation. 


Hancock’s predictor corrector scheme. 


Consider the simple linear one dimensional convection equation 

&- + a^-=0 (4.1) 

3t 5x 

We solve this numerically by a MUSCL scheme in a network t=nx , x = jh where integral 
values of j correspond to cells and half integral values to cell boundaries or faces. 


According to the original MUSCL scheme (van Leer (1979)) the following steps are exe- 
cuted to update the cell interface values from time t=nx to t=(n+l)x 

a) ptjj/, = p 9 + — Ap " . Space extrapolation to cell interfaces 


b) pp/{ 2 = PjJ.i/ 2 + — Time extrapolation to t = (n + Vi) x 


where p£i /2 = a 


Apj 

Ax., 


C ) P^'/2 1 = P jh'/2 + X P jh'/2 

Face value of p at t = (n+l)x to be used in interface differencing. 


d) pf 1 = pj - a^ 

Control volume updating of cell mean value. 


n+'/2 ^ n+'/2 

P jh 1 / 2 P y-Vi 


e) Apj n+1 = ppl - pp} 2 New slope 


rw-1 


f) Limit Ap/ 
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The iterative process needed to carry out steps a), b) & c) is simplified by Hancock by the 
following two step method. 

a) p ? +I = p? - a — Ap" Extrapolation of cell mean values forward in time, using PDE. 


n+'/2 L 


n+72 _ 

b) Pi+'/ 2 - 


Pj + Pj 


n+1 


+ y A p? Select upstream cell values at n +'/2 


c) p£J^ = p" +1 + y Ap" Select upstream face values at n+1 


Steps d) e) and f) follow as before. 


This simplified procedure is extended to the one dimensional Euler equations by van Leer 
in van Albada, van Leer and Roberts (1982). 


The van Leer version of the higher order Godunov method comprises two main steps, 
reconstruction and evolution. In the first step the data. obtained from the most recent difference 
calculation are fitted by piecewise linear representation. The fitting is straightforward over 
regions in which the data are smooth and monotonic. At discontinuities, or local extrema the 
procedure must be modified to preserve monotonic behavior in further computation. At such 
locations constraints are imposed and the details of these will be given. 


The second evolution step uses an upwind procedure developed from the MUSCL scheme. 


Let W represent a complete set of state variables such as p, pu, p or p, u, p. 
The initial distribution of W is represented in a piecewise linear fashion by 


(AW)? 

W n Ol) = W? + (Tl — Tlj )-— 3 - 


4.2 


where 
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(Au)j n = c ■ ave 


n n _ n n n_ n n 
9j Qj-i 


4.3 


and 


(Ap)j n = pj n - ave 


2 (Pz j+ i-P f) ^ 2 (Pj B -P^i) 


(Pzj+i + Pj) (Pj+P*l) 


4.4 


where ave (a, b) is an average specified in van Albada et al (1982). 


Then 


aw 


dn 


" (AW)j n 
i At| 


4.5 


and 


aw 


at 


can then be found from the Euler equation. 


The cell averages can now be determined, at time t n+ y 2 leading to boundary values. 


wp* = W? + — 
— 3 — 3 2 


\\t n+'A _ \\7 n ± */2 
— ( j±‘/2); ~ Wj 


aw 


at 

AW 


4.6 


4.7 


The vector U of conserved variables ( p,pu, e) can then be found from 


U (j+ ./:); = U 


w 


n+!/2 

(j±‘/2)l 


4.8 


These in tum give the values of time centered fluxes at cell boundaries. 


5. Huynh’s Improvement of the MUSCL Scheme. 


The generalization of Hancock’s improvement of the MUSCL scheme by Huynh (1985) is 
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carried out in three stages, discretization, slope constraint and Riemann solver. In the discretiza- 
tion procedure the Eulerian equations for one dimensional unsteady flow, in conservation form, 
are written as difference equations using a mid-point rule. Here, the fluxes at adjacent cell boun- 
daries need to be evaluated at a half time step, after breakdown, by solving a Riemann problem. 
The required initial values of the conservation variables on either side of a cell boundary are cal- 
culated from linear representations in space and time. The Riemann problem uses Roe flux- 
difference splitting (Roe (1981) and Roe (1986)). 

Slope constraint is needed to maintain the monotonic character of the numerical solution 
for the conservation variables, as this is advanced in time. The increased accuracy of the 
Godunov method, consequent on its upgrading to second order, is achieved at the the price of 
generation of local maxima or other violations of monotonic behavior within the finite difference 
network. The slope constraint stage describes the procedure for correcting these violations. 
Huynh’s description of this in algebraic form is complicated and details can be found in his 
paper. It is somewhat more general than the procedure shown graphically by van Leer (1979) 
which is reproduced later. 

The approximate Riemann Solver proposed by Roe (1981, 1986) is the principal tool used 
by Huynh in the third stage of his scheme and is presented here in Roe’s compact form. 

In essence the improvement in the MUSCL scheme by Huynh rests on a recognition by 
Hancock that, in the two half-time steps taken in MUSCL, the Riemann problem only needs to 
be solved at one of these steps. In Huynh’s formulation this is done at the half time step. 

Discretization 

Huynh uses the Eulerian equations of one dimensional unsteady flow in both conservation 
and primitive form. For a perfect gas, with constant specific heat ratio y the conservation equa- 


tions are 
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9U 9F(U) 

at ax 


(5.1) 




r 

p 


r *\ 

pu 

where 

U = 

pu 

e 

L J 

F = 

pu+p 

>+Pl 


with 


e = p/(y- 1) + j pu 2 

The primitive variables are 


(5.2) 


V = (p, u,p) T (5.3) 

From the governing equations for the primitive or conservation 

variables we can derive the equations in characteristic form (see 1.3 in Holt (1984)). 


To discretize Eq. (5.1) introduce a uniform mesh defined by Xj = jh, j = 0, 1, 2 • • • 
t n = n x n= 0, 1 , 2 • • • where h is the mesh size and x the time step. The cell boundaries are at 

x i = (j + — ) h. Let Uj n be an approximation to the average value of U in the jth cell at time t n 
*2 2 


U; n = — f U (x, t n ) dx 

J h J 


l 

X)+ 2 


(5.4) 


Given U" we wish to determine Uj l+1 at the end of time step x, subject to the CFL condition 


max ( I Uj I + Cj) 


L J 


-< 1 
h 


(5.5) 


where Cj is the local speed of sound. 
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Differencing Eq (5.1) using the mid-point rule gives 


U; n+1 = Uj" + i- 

J J h 


n+ 


1 


n+ 


F { -F ± 
J 2 J+ 2 


(5.6) 


n+ . 


To solve Eq. (5.1) we need to determine the fluxes F” , r at the half time step from the initial 

J+ 2 

data Uj n or, equivalently, from V". 

In each j th cell we represent the primitive variable V(x, t) by a linear function 


Rj(x,t) = Vj + (x - xj) Sj + (t - t n ) Tj 
where Sj = —■ (xj, t n ), Tj = (xj, t n ) 


(5.7) 


dx V ' T * ” ‘ J 9t 

Sj is known from the initial distribution of V and Tj can be determined from the primitive form 
of Eq (5.1) 


V-(A p)jSj 

where 


up 0] 



0 u 1/p 

,0 TP u 


At the interface j + Vi Eq (5.7) gives 


(5.8) 


(5.9) 


Rj (x ,,t" +2 ) = V j + |s j + iT j (5.10) 

S+2 2 2 

Eq (5.7) gives values of the primitive variables on the left side of the interface j + V 2 , resulting 
from analysis in cell j. Similar analysis in cell j + 1 gives the values of the primitive variables on 
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the right side of the interface j + Vi. These calculations provide the initial data for a Riemann 
problem to be solved on interface j + Vi at time (n + Vi) x. To this end Huynh uses Roe’s 
Riemann solver. 

Riemann Solver 


Huynh uses a generalized version of Roe’s method for solving Riemann problems which 
arise in Godunov’s method (first or second order versions), approximately. We present a sum- 
mary of Roe’s version (1981), Roe (1986). 

Using Huynh’s notation Roe uses a locally linearized form of Eqs (5.1) 

U t + AU X = 0 (5.11) 

where A is a constant Jacobian matrix 8F/3U. If U L and U R are the interface values of U on 
either side of (j + Vi)Y\ the flux difference across the interface can be written 

F R -F L = Za k ?l k e k (5.12) 

where ( e k ) are the right eigenvectors of A. Each term on the right of Eq. (5.12) gives the effect 
of the k th wave crossing the interface, with a k denoting its strength and X k ( eigenvalue of A) its 
speed. The flux at the interface (j + l /z)h can be computed either from summation over negative 
or over positive wave speeds. Roe takes the average of the two evaluations to give 

F ^ (U L , U R ) = Vi (F l + F r ) - E ak I A, k I e k (5.13) 

*2 

To apply (5.13) to the present non-linear problem we use a matrice A (U 2 , U R ) the eigenvalues 
and eigenvectors of which not only satisfy (5.12) identically but also 

^R _ ^L = ^ a k e k (5-14) 

k 

Under these conditions the method, for non-linear problems, gives the exact solution when U L 
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and U R are on opposite sides of a shock wave or contact discontinuity. 
The expressions for a k , A. k and e k , derived by Roe (1981) are 



f ^ 

1 


1 


1 

*1 = 

u - a 
h - ua 

w J 

. «2 = 

U 

Vi u 2 

k. j 

. «3 = 

u + a 
h + ua 


X, = u - a, X 2 = u , A, 3 = u + a. 


1 1 2 

aj = — - [Ap - paAu] , a 2 = -r[a Ap - Ap] , 

2£r a 


a 3 = — T - [Ap + paAu], 
2a 2 


where 


(5.15) 


(5.16) 

(5.17) 


p 2 - Pl Pr ’ 


1A |A 

Pl Ul + Pr u r 
Pl + Pr 


Pl 2 h L + Pr^r 
Pl 2 + Pr 


VH 


h = 


(5.18) 


a 2 = (Y-l)[h-y u 2 ]. 

When we introduce (5.15) thro, (5.18) into (5.13) we obtain an explicit expression for the 

n+ n+ \- 

required fluxes F { and F j in 5.6. 

J_ 2 J+ 2 

Monotonicitv constraints 
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On completion of a time step in the MUSCL scheme (original or improved version) the 
representation of the unknowns by linear functions in successive cells must be checked to ensure 
that no local violations of monotonic behavior take place. This condition is most simply 
expressed by van Leer (1977). Namely, in terms of Eq. (5.7) at time t = t n , the linear x function 
must not take values outside the range spanned by the neighboring mesh averages. 

Fig. 1 ( taken from van Leer (1977)) shows the three possible violations. In section (1) of 
this a local decrease at 0 violates monotonicity. To remedy this the slopes in (-1, 0) and in (1, 2) 
cells are reduced to zero while the slope in (0,1) is reduced as shown by the solid line. In section 
(2), where a local extremum is attained at boundary 1 the slope in (0,1) is reduced to zero. In 
section (3) the sign of the slope in (0,1) is opposite to the neighboring slope signs in (-1, 0) and 
(1,2) and is therefore reduced to zero. These corrections can all be expressed in algebraic form 
and are detailed in Huynh (1995). 



Fig. 1 


6. TVD and ENO Schemes 


A later and somewhat different study of higher order Godunov methods was initiated by 
Harten (1983). This started with introduction of the Total Variation Diminishing (TVD) 
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schemes and followed by this broader Essentially Non-Oscillatory (ENO) schemes. 

Davis (1984) provides an outline of the TVD scheme applied to the one dimensional wave 
equation. 

Consider the Initial Value Problem 


u t + f(u) x = u t + a(u) u x = 0 


a(u) = 


df(u) 

du 


( 6 . 1 ) 


U(X, 0) = u 0 (x) - co < X < °° 


where u 0 (x) is of total bounded variation. 


We seek a weak solution to this problem with the following properties. 


(1) No new maxima or minima in u(x) may develop 


(2) The value of any local minimum does not fall while that of any local maximum does not 
rise. The Total Variation of the solution is defined by 


TVu(x.t) = sup £ I u(x k+1 ,t) - u(x k ,t) I 
k 

where the supremum is taken over all partitions of the real line. 

The monotonicity property requires that the total variation in x of u(x, t) does not increase 
with t so that 


TV{u(t 2 )} <TV{u(t,)} for all t 2 >t,. 


Now consider an explicit finite difference scheme in conservation form applied to IVP - 6.1 
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U n+1 = L • U n (6.2) 

where L is an explicit finite difference operation. 

The scheme is said to be Total Variation Diminishing if 

TV (U n+1 ) =TV{L • U n } <TV {U n } (6.3) 

Furthermore, a scheme is Monotonicity Preserving if L has the following property. If U n is 
a monotonic mesh function then L • U n is also monotonic. Finally Harten proved theorem: 

A Total Variation Diminishing Scheme is Monotonicity Preserving. 

Harten (1986) notes that TVD schemes, which include the MUSCL scheme, are second 
order accurate only in the L] sense. To increase this accuracy he introduces the Essentially Non- 
Oscillatory (ENO) in which the Total Variation requirement on the unknown, at the end of a 
time step, is relaxed in the reconstruction phase. ENO schemes can be accurate to any finite 
order. 

Conclusions 

Extensions to second-order accuracy of Godunov’s method for solving Gas Dynamics 
equations numerically are reviewed. Particular attention is given to the MUSCL scheme and 
its recent improvement by Huynh due to a simplification first noticed by Hancock. 

The author is indebted to Colella, Hancock, Huynh, Roe and van Leer for providing 
original papers and to van Leer for enlightening discussion. The author thanks M.Y. Hussaini 
for arranging visits to ICASE, where most of the paper was assembled. 
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